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Abstract 

Quarkonium at finite temperature is described as an open quantum system whose dynamics are 
determined by a potential Vr(x) and drag coefficient rj, using a path integral with a non-local term. 
Path-integral Monte Carlo calculations determine the Euclidean Green function for this system to 
an accuracy greater than one part in a thousand and the maximum entropy method is used to 
determine the spectral function; challenges facing any kind of deconvolution are discussed in detail 
with the aim of developing intuition for when deconvolution is possible. Significant changes to 
the quarkonium spectral function in the 15 channel are found, suggesting that any description 
of quarkonium at finite temperature, using a potential, must also carefully consider the effect of 
dissipation. 
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I. INTRODUCTION 


At both zero temperature and above deconfinement, quarkonium is an ideal probe of 
QCD. At zero temperature its mass M 3> Aqcd suggests the effective field theory approach 
of NRQCD [1], where the heavy quark mass is integrated out. In effect, the heavy quark mass 
sets the momentum scale of QCD for states using this description, allowing perturbative re¬ 
sults to work well. The situation becomes more complicated at finite temperatures. NRQCD 
can be extended to finite temperature and used to examine the break-up of quarkonium [2] . 

An alternative point of entry is to treat quarkonium above deconfinement as an open 
quantum system (OQS)|l]. This departs from any attempt to describe quarkonium with 
perturbation theory, and the parameters describing quarkonium have to be determined sep¬ 
arately, either in calculations based on first principles or from experimental measurements. 
However, this description has the strength of being independent of some of the scale hier¬ 
archies. One would suspect that the OQS approach is at its most useful for charmonium 
and highly excited states of bottomonium, where the binding energy is small and the con¬ 
fining term of the Cornell potential is important for describing the spectrum of states at 
zero temperature. Perhaps the greatest strength of the OQS approach is that the approach 
to thermalized yields of quarkonium is completely natural, which is not necessarily the case 
for some potential models. Since the first, simplest descriptions of quarkonium in this way, 
the treatment of heavy quarks at finite temperature has developed significantly, thanks the 
work of several authors ISHT]. 

This paper describes the numerical determination of the spectral function for an open 
quantum system, which has the potential term and the diffusion coefficient matched to 
quarkonium at the temperatures reached in heavy-ion collisions. The procedures used can 
be generalized for different densities of states for the thermal bath, for different couplings, 
and for non-trivial correlations of the bath’s force on the heavy quark and the anti-quark. 
In Section the path-integral Monte Carlo algorithm for this open quantum system is 


described. In Section intuition is developed for when phenomenologically significant de- 
convolution of spectral functions is possible, with this having some implications for results 
using lattice QCD calculations; and finally, the maximum entropy method is used to de¬ 
convolve the results from Section into quarkonium spectral functions, and the non-trivial 
relationship between diffusion and the destruction of the J/ip state can be determined. In 
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this way, the often competing effects of drag and momentum diffusion on J/ijj survival rates 
can be simultaneously considered, as was done in Pi- 

While this procedure may seem cumbersome, it is in fact relatively cheap computationally. 
The most costly step is the calculation of the correlation functions to sufficient accuracy so 
that deconvolution yields signihcant results; work on high-order estimates of the action in 
Equation can speed this computation signihcantly as it has for simple actions. 


II. EUCLIDEAN CURRENT-CURRENT CORRELATORS AND DISSIPATIVE 
EFFECTS 


Much of the groundwork towards examining open quantum systems with path integrals 
was already done by Feynman and Hibbs when they considered influence functionals HDl; 
Caldeira and Leggett used path integrals to describe a system with a specihc thermal bath 
[TT] while Grabert, Schramm, and Ingold generalized Caldeira and Leggett’s results [I2]. In 
|1] , these techniques were applied to imaginary-time Green functions. To review: one starts 
with the action for a heavy degree of freedom with position x interacting with a light degree 
of freedom described by R. This light degree of freedom is often a collective mode of a gas 
or of condensed matter, and not just a value for position. The action for this system, using 
a harmonic approximation and minimal coupling, is 

+ V{x) + + ^moj'^R^ — CxR 

gives the influence functional 



{xf,T\xi,0)^^^ = j exp 

G2 


dr 


-I- V{x) 


2ma;sinh(a;r) 


x(r ) cosh(c(;(r — T )) / ds x(s) cosh(c(;s) 


( 1 ) 


after Ending {xf, Rf\xi, Ri) and integrating over Ri and Rj. When there are mul¬ 
tiple light degrees of freedom, the final term in the integral becomes the sum 
eosh(a;j(r — r')) JJ" ds x{s) cosh(a;js). The system becomes dissipative 
in the limit of an infinite number of light degrees of freedom of masses rrii and frequencies 
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oji- the density of states 


C‘^{u)pd{u:) = 


( 2 ) 


if a; < 

TT 

0 if w > f2 

in the limit f2 —)■ oo was examined in m and fonnd to lead to a path integral which yields 
the Langevin eqnation in the classical, high-temperatnre limit. Performing the integral over 
these states in Eqnation [T] leads to the integral 


du 


du r dv - ^)) 


sinh(a;r) 


n 


= — du{x{u)y — 


du 


Iq C(;sinh(a;r) Jq 


du / dn i;(M)i;(n) sinh(a;(r — m)) sinh(a;n), 


which forces a renormalization of V{x), and gives the imaginary-time Green fnnction for 
this system: 


Gred{Xf,Xi,T) = Vx exp i - 


du 


^Mxiuf + Vr{x{u)) - ^ 


dv x{u)x{v) log 


sin(f^; 

sin(f^^; 


Note that only the potential is renormalized, and that the hnal term introdnced by the 
integration over the bath of particles is translationally invariant, meaning that no hnite, 
x-dependent term has been added to the path integral. This imaginary-time Green fnnction 
can be made periodic in r with period (3 nsing the method of images to make a qnantity 
related to results from hnite-temperature lattice calculations. 

For a general potential, numerical methods for determining this path integral must be 
developed. Path-integral Monte Garlo techniques have been developed for condensed matter 
systems; Geperley reviewed these methods as they concern liquid helium [13]. For quarko- 
nium, a relatively simple algorithm is used: the path is discretized to 2^^ -|- 1 equally-spaced 
times, and free-particle paths are sampled using the bisection method. These paths are 
reweighted according to the expression above. The double integral in the exponential poses 
some problems because it is improper; both the “square” and “triangular” regions in the 
double integral must have their measures determined analytically so that sufficient precision 
may be achieved. 
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The results for G{t) have been calculated for the Cornell-like potential 


l.bois 


V{r) = < 


-f O.Scrrmin + (7^ + O.Scrrmin)^ 

- iiiixx ' min ' min 

if r < rmin, 

-as/r + ar lir > rmm 


( 3 ) 


with Tmin chosen to be small, 0.4 GeV“^ (here, all units are in GeV), as = 0.499, and 
cr = 0.16 GeV^. The purpose of using this piecewise function is to simplify dealing with 
the divergence of the Gornell potential as r —0. The drag coefficients rj =0, 0.0729 GeV, 
and 0.1458 GeV were used, corresponding to the spatial diffusion coefficients 2ttDc = 00 , 
5, and 2.5, respectively, at T = 285 MeV. The results for these values were shown in |1]; 
here, in Figure [T| we show results for the same values but for a large range in r, which is 
necessary for extracting the spectral function. We should also note that some degrees of 
freedom of the heavy quark are not represented in this path integral, namely, spin and color. 
The potential term is therefore the thermal average of this quantity over spin and color 
states. In [6], Lindblad equations describing the evolution of color singlet and octet states 
are determined for heavy quarkonium, starting from perturbative QGD for heavy quarks. 
Our path integral, on the other hand, requires color-averaged potentials as an input for the 
path integral Monte Garlo calculations. 


III. THE QUARKONIUM SPECTRAL FUNCTION AND THE CHALLENGES 
FACING DECONVOLUTION 


The imaginary-time Green function is related to the spectral function p = —Im{GR} 
through a Laplace transform: 


G{t) = / exp(—ci;r)p(ci;)(ia;; 


( 4 ) 


similarly, the imaginary-time hnite-temperature Green function is this Laplace transform 
made periodic through the method of images: 


f cosh(n;(r-/?/2)) 

G{r,^)= / - . , . a/ON - p[^)duj. 


( 5 ) 


sinh(n;/?/2) 

Gomplex analysis helps here: the Fourier transforms of the various Green functions are con¬ 
veniently related to each other; as a function of complex cu, causality translates to analyticity 
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FIG. 1 (Color Online): From [3], G{t) with and without dissipative effects. 

of Gn in the upper-half plane; and the relation between Gr and Ga in o; is a generalization 
of the fluctuation-dissipation theorem [T3] . 

Unfortunately, when working with numerical results and not analytic expressions for G{t), 
these results are not helpful. The reason for this is simple: signihcantly different spectral 
functions often differ by very little in their Laplace transforms. This is not obvious to 
theorists who work with analytic results so it is necessary to illustrate this with an example: 
consider two spectral functions, pi and p 2 , plotted in Figure]^ The widths differ by a factor 
of hve, from 0.01 GeV for pi to 0.05 GeV for p 2 - This changes the lifetime of the state from 
~ 20 fm/c to ~ 4 fm/c. The lifetime of the state represented by pi is long compared with 
the timescales of the heavy-ion collisions at RHIG and the LHG; it represents a state whose 
yields would be largely unaffected, while the state represented by p 2 would be signihcantly 
suppressed. 

Lattice QGD calculations of quarkonium correlation functions have determined G(r, (3) 
for temperatures between = 175 MeV and 2Tc, by determining the correlation functions 
of composite operators 

G{r,^) = J , ( 6 ) 

where and T^ = 7 ^, 747 ^^. This function corresponds to the transform of the 
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FIG. 2 (Color Online): The spectral functions pi (with width cr = 0.01 GeV) and p 2 (with width 
a = 0.05 GeV). 

spectral function shown in Equation]^ For f3 = (l.ST^)”^, this transform of pi and p 2 can 
be computed; Figure shows the relative difference AG/G = (G 2 — Gi)/G 2 for values of r 
from zero to 0.5/?. 

The relative difference has a maximum of about one part in a thousand. If the full range 
for r were equally sensitive to changes in this bound state, this would suggest that three 
results for G(r,/?), equally spaced on this range, would need an accuracy of a few parts in 
ten thousand to determine which spectral function hts the results best. Constraining the 
spectral function in the vector channel with sum rules is not likely to help much because 
they provide only a few constraints to a continuous function. 

Lattice QCD calculations at specihc temperatures are performed on periodic lattices and 
are used to determine hnite-temperature Green functions as in Equation not However, 
if another calculation can determine imaginary-time Green functions as in Equation at 
large r, it will be sensitive mostly on the shape of the ground state and the resulting Gi and 
G 2 should differ by a greater relative AG. We take advantage of this fact by calculating this 
Green function for a dissipative system over the range shown in Figure]^ Here, |AG|/G has 
a maximum of about 1%, suggesting that numerical deconvolution will be less daunting. 
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FIG. 3 (Color Online): The relative difference AG/G = (G 2 (t, /3) — Gi{t, /3))/G 2 {t, /3) between 
the Green functions obtained from pi and p 2 - 



FIG. 4 (Color Online): Gi(r), G 2 {t), and |AG|/G, plotted together. 







If G'(r) is determined with sufficient accuracy, deconvolution can be attempted. Again, 
the inverse Laplace transform is analytically easy but numerically non-trivial. Working with 
tens of data points for G(r), no matter the accuracy, underdetermine any reasonable dis¬ 
cretization of the spectral function. Consider y^-minimization in this situation: the spectral 
function is discretized into ~ 1000 points and ~ 10 data points are htted, unnaturally small 
values of can be achieved. The result for p{uS) usually appears choppy in this situation. 
The principle of maximum entropy suggests that the best £t for piuj) is not the £t with 
the smallest but the £t with a reasonable value for but is also constrained by the 
information entropy 

^ = log (Pi/- (Pi - ’ ('") 

i 

where pi are the discretized values for p{u:) and cTj are the priors for pi. This constraint on 
the x^-minimization is often cast in terms of Bayesian principles, but it can also be thought 
of as an assumption of smoothness for the deconvolved function, far less controversial in 
physical situations. The sum of cind I dehne the energy function to be minimized: 

E{pi) = X^{pi) + Oil. ( 8 ) 

The coefficient a determines the relative importance of x^ and I in the minimization; it is 
often chosen to make x^ roughly equal the number of data points. Gallicchio and Berne 
point out that the best value for a can also be determined with Bayesian logic na. Our 
choices for a will be based simply on whether or not they yield reasonable values for 

A suitable algorithm for the minimization of this function is simulated annealing: 1.) 
random steps are considered in the space of values for p*, , changing the value of E from 
Ei to Ef, 2.) any step that decreases E is taken while any step increasing E is taken with 
probability P = exp{—{Ef — Ei)/T), T being a temperature chosen to be large initially, 
and 3.) the process is repeated with T lowered until a minimum temperature T is reached. 
When a function is multi-modal as is possibly the case, simulated annealing is far more 
likely to be successful at finding absolute minima than the biconjugate gradient method, or 
any other slope-following method. 

At this point, a test of this method of deconvolution would be useful. We can test this by 
going in the opposite direction: starting with a given spectral function, hnding its Laplace 
transform (with random Gaussian error added), and deconvolving with the maximum en- 
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a=0.01 

a=0.1 

a=l 

1 st peak 

17.15% 

14.84% 

3.58% 

2 nd peak 

22.85% 

10.29% 

4.65% 


TABLE I: The discrepency between the true spectral function and the results of deconvolution 
shown in Fig. 

tropy method. The spectral function 

_(u-0.33)2 _(i.j+0.33)^ 

(j(0J^ = 1.2(e 2*0.0174^ — e 2*0.0174^ ) 

2 2 (9) 

+ 0.2(e 2 * 0 . 108 ^ — e 2 * 0 . 108 ^ j 

mimics the form that we expect for the quarkonium spectral function at finite temperature. 
The data set 

G{Ti) = j a{u)e~‘^'^du + AG* (10) 

contains 15 points ranging from r = 0.5 GeV“^ and r = 9.6 GeV“^, and AG* is random 
Ganssian noise with standard deviation 10“® GeV. In Fignrej^ the resnlts of the maximnm 
entropy method with different valnes of a are compared to the original spectral fnnction. 
In Table 1, the error in the widths of the peaks in the spectral fnnction are shown for the 
different valnes of a, showing agreement within a few percent when a = 1. 

The resnlts from applying the maximnm entropy method to the resnlts for the qnarkoninm 
Green fnnctions are shown in Fignrej^ It has resolved what shonld be a Dirac delta fnnction 
in the spectral fnnction withont dissipation down to a width of 0.25 GeV. The spatial 
diffnsion coefficient 2nTD = 5 corresponds with an increase of the width to 1 GeV, signifying 
the state having a lifetime of 0.197 fm/c. The resnlts of the simulated annealing for 27rTD=cx) 
with different a valnes can be seen in Fignrej^ It is observed that the resnlts are not sensitive 
to the valne of a over a considerable range. 

We end this section by noting the physics behind the changes of the spectral fnnction in 
Figure First, the centroid of the gronnd state peak is shifted to high uj with increasing 
rj. This snggests that the mean radins of the gronnd state increases, nnder the inflnence of 
momentnm transfer from the medinm. This conld honestly have been seen rather qnickly, 
by examining the slope of G(r) in Fignre A far less trivial resnlt from the maximnm 
entropy method is the increasing width of the gronnd state peak, indicating a decrease of 
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FIG. 5 (Color Online): The spectral functions deconvolved from the test data by using different 
values of a. The corresponding values of are shown in the legend. 



FIG. 6 (Color Online): The spectral functions deconvolved from the results in Figure The 
corresponding values of x^ are shown in the legend. 
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FIG. 7 (Color Online): The spectral functions deconvolved from the results where r/ = 0 in Figure 
[^with different values of a. The corresponding values of are shown in the legend. 

the lifetime of the ground state. 


IV. SUMMARY 

Deconvolution of spectral functions was examined with an emphasis on the building of 
intuition for quarkonium near T^. Results from a treatment of quarkonium above deconfine¬ 
ment as an open quantum system yielded a set of spectral functions showing the effect of 
increasing rj, and these results were found to be robust for a range of values in a. 

While this paper dug very deeply into the issues related to deconvolution, more work 
is required. In particular, only flat priors were used here; the role of the prior must be 
examined carefully. Finally, temperature-dependent potentials will be used in future work, 
making possible strong statements about quarkonium spectral functions including multiple 
temperature-dependent effects. 
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